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Abstract 


Here we demonstrate that the activity of neural ensembles can be 
quantitatively modeled. We first show that an ensemble dynamical model 
(EDM) accurately approximates the distribution of voltages and aver¬ 
age firing rate per neuron of a population of simulated integrate-and-fire 
neurons. EDMs are liigh-dimensional nonlinear dynamical models. To 
faciliate the estimation of their parameters we present a dimensionality 
reduction method and study its performance with simulated data. We 
then introduce and evaluate a maximum-likelihood method to estimate 
connectivity parameters in networks of EDMS. Finally, we show that this 
model an methods accurately approximate the high-gamma power evoked 
by pure tones in the auditory cortex of rodents. Overall, this article 
demonstrates that quantitatively modeling brain activity at the ensemble 
level is indeed possible, and opens the way to understanding the compu¬ 
tations performed by neural ensembles, which could revolutionarize our 
understanding of brain function. 
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1 Introduction 


If we observe a fluid at the molecular level we see random motions, but 
if we look at it macroscopically we may see a smooth flow. An intriguing 
possibility is that by analyzing brain activity at a macroscopic level, i.e., 
at the level of neural ensembles, we may discover patterns not apparent 
at the single-neuron level, that are as useful as velocity or temperature 
are to understand, and predict, the motion of fluids. 

Technology frequently drives science. For instance, thanks to the de¬ 
velopment of microelectordes in the 1930’s, we now know with exquiste 
detail computations performed by single neurons. We are now experienc¬ 
ing a dramatic increase in our capacity to monitor the activity of larger 
and larger populations of neurons with higher and higher spatial and 
temporal resolution. These new ensemble recordings may soon allow us 
to uncover crucial computations performed by neural ensembles. 

Here we present results of developing and evaluating a mathematical 
model and estimation methods to characterize the activity of ensembles 
of neurons from electrophysiological data. 

Section[2]reports the evaluation of an ensemble dyanmical model (EDMs). 
And EDMs is a high-dimensional nonlinear dynamical models. To esti¬ 
mate its parameters it is convenient to reduce the number of parameters. 
We escribe a dimensionality reduction method for EDMs in Section [3] 
Section [4] presents a maximum-likelihood method to estimate parameters 
of EDMs. Finally, in Section 0 we show that EDMs can accurately ap¬ 
proximate high-gamma electroencephalographic activity evoked by pure 
tones in the auditory cortex of rodents. 
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2 Building EDMs 

We wanted to learn how to build Ensemble Density Models (EDMs), dy¬ 
namical models of the state variables (e.g., trans-membrane potential, 
time since last spike, etc.) of a population of identical neurons, starting 
from a dynamical model of the state variables of a single neuron. For In¬ 
tegrate and Fire (IF) model of single neurons, an EDM should provide the 
ensemble probability density function (pdf) p(v,t), from which to com¬ 
pute the probability of finding a neuron in the ensemble with with a given 
trans-membrane voltage v at time t (i.e., P(v,t) = p(v,t)dt). It should 
also provide the average firing rate per neuron in the pop ulation, r(t). To 
constr uct EDMs we chose the methodology described in lOmurtag et alj 

Em. 

To evaluate the EDM we compared its outputs (p(v,t) and r(f)) with 
those derived from the direct simulation of a population of 9,000 IF neu¬ 
rons. The value of the density function p(v,t) derived from the direct 
simulation was the proportion of IF neurons having a voltage v at time 
t, and the average firing rate per neuron r(t) derived from the the direct 
simulation was the proportion of cells in the population at time t with 
voltage at threshold. 

An EDM is driven by an external excitatory current and modulated 
by an external inhibitory current. In addition, every cell in the EDM 
receives inputs from G other neurons in the population. A fraction / 
of these G inputs is excitatory, and the remainder are inhibitory. These 
intra-population inputs act as feedback mechanisms to the EDM. The 
mathematical representation of an EDM for a population of IF neurons 
is given in Eq u ation s [l][TTJ modified from Equations 26, 39 and 47 in 
lOmurtag et al.l i20C)dl . External excitatory and inhibitory currents ap¬ 
pears as CTe(t) and o-°(t), respectively, in Equation [3l and excitatory and 
inhibitory feedback are given by the terms Gfr(t) and G( 1 — f)r(t), re¬ 
spectively, in Equation [3] 


dp 

dt 


(v,t) 

r(t) 




dJ . , 

J(v = 1, t) 

-7 vp(v,t) 

+[cr°(f) + Gfr(t)] f p(v',t)dv' 

J V- 


-K°(f) + G(l-/)r(f)] 


1 , 




p{v', t)dv' 


(1) 

( 2 ) 


(3) 


We first verified that for a single population the outputs of the EDM 
matched those of the direct simulations (Section 12.111 . We next built a 
network of excitatory and inhibitory populations, and again compared 
the outputs of the EDM and those of the direct stimulation ('Section l2.2l) . 
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2.1 One Population 

2.1.1 One Population of Independent Neurons 

The top panel in Figure [T] shows the ensemble probability density func¬ 
tion, p(v, t ), calculated by integrating the differential equation of an EDM, 
Equation [l] The bottom panel shows and approximation to this pdf ob¬ 
tained from the histogram of voltages of a direct simulation of a population 
of 9,000 IF neurons. Note the large similarity of the pdfs at all time points. 

The black line in Figure [2] shows the average firing rate per neuron, 
r(t), calculated using the EDM, Equation [5] The grey line shows the 
average firing rate per neuron calculated from a direct simulation of a 
population of 9,000 IF neurons, as the proportion of cells with voltage at 
threshold. Note the almost perfect match between these firing rates. 

FigureEfs as FigureEJbut for a step input current that jumps from 0 to 
800 impulses per second at time 1 = 0. Note that by 0.4 seconds after the 
step in the input current the average firing rate per neuron has reached 
a new steady state around 10 impulses per second, as revealed by the 
EDM (black line in Figure[3| and by the direct simulation of a population 
of 9,000 IF neurons (grey line, in Figure 0). Figure U shows the pdf 
of voltages at this new steady state calculated by the EDM, Equation E] 
(black line in Figure [4]) and approximated using the histogram of voltages 
at 0.4 ms of a direct simulation of 9,000 IF neurons (grey line in FigureE]). 

2.1.2 One Population with Feedback 

The previous Figures showed results from the simulation of a single pop¬ 
ulation of independent neurons. Figure [5] is as Figure [2j but shows the 
average firing rate per neuron in a population where each neuron receives 
excitatory inputs from ten other neurons in the population (i.e., G=10 and 
f=l in Equation El- For comparison, the dashed line shows the average 
firing rate per neuron from the population without feedback. 

Figure E] demonstrates the effect of inhibitory feedback in the popu¬ 
lation of Figure [5] by changing the fraction of inhibitory input neurons 
from 0% to 80% (by changing f=l-0 to f= 1-0.8, and still using G=10 in 
Equation El) ■ 

2.2 A Network of Populations 

The previous sections evaluated EDMs in a single population of neurons. 
Here we evaluate EDMs of excitatory and inhibitory populations combined 
in the network of populations shown in Figure [7] The network is driven 
by an excitatory input to the excitatory population. This population has 
excitatory feedback (i.e., G=5 and f=l in Equation El). and its average 
firing rate per neuron output, scaled by a constant W e i = 50, drives the 
inhibitory population. The inhibitory population has inhibitory feedback 
(i.e., G=5 and f=0 in Equation El, and its average firing rate per neuron 
output , scaled by a constant Wi e = 15, modulates the activity of the 
excitatory population. 

Figures El and E] show the activity in the excitatory and inhibitory 
populations, respectively in the network of Figure [7] The upper panels in 
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Figure 1: Ensemble pdf, p(v,t), for a population of IF neurons in response 
to a sinusoidal input, calculated by integrating the differential equation of an 
EDM (Equation [l] top panel) and approximated by direct simulation of a pop¬ 
ulation of 9,000 IF neurons (bottom panel). Neurons in this populations were 
independent of each other (i.e., they had no feedback; G=0 in Equation [3]) . 
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Figure 2: Firing rate of a population of neurons, in response to a sinusoidal 
input, calculated by direct simulation of a population of 9,000 IF neurons (grey 
line) and by integrating a population equation (black line, Equation[3]). Neurons 
in this populations were independent of each other (i.e., they had no feedback; 
G=0 in Equation [3]). 





























Figure 3: Firing rate of a population of 9,000 IF neurons, in response to a step 
input at time zero. Same format as in Figure [2] 
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Figure 4: Ensemble pdf, p(v,t) at 0.4 seconds in the neurons of the population 
of Figure [3] in response of a step input at time zero. This pdf was calculated 
by integrating the differential equation of an EDM (Equation |T| black line) and 
approximated using the histogram of the voltages of the simulated neurons (grey 
line). 
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Figure 5: Firing rate of a population of 9,000 neurons, in response to a sinusoidal 
input, as in Figure [2j but each neuron in this population is connected with ten 
presynaptic neurons, and all of these neurons are excitatory (i.e., G=10 and 
f=1.0 in Equation [3]). 
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Figure 6: Same as Figure [5] but for a population of neurons where 80% of the 
presynaptic neurons to a given neuron are inhibitory (i.e., G=10 and f=l-0.8 in 
Equation [3j. For comparison, the dashed line shows the average firing rate per 
neuron in the population with feedback but not inhibition of Figure 0] 
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Figure 7: Simulated network of two populations of IF neurons. A sinusoidal 
input, er°(f), was applied to the excitatory population. This population con¬ 
tained excitatory feedback (each neuron in the population received excitatory 
input from five other neurons in the population; G=5, f=l in Equation [3jl. 
The average firing rate per neuron of the excitatory population, scaled by the 
coefficient W e i = 50, was the excitatory input for the inhibitory population. 
This population contained inhibitory feedback (each neuron in the population 
received inhibitory input from five other neurons in the population; G=5, f=0 
in Equation[3]). The average firing rate per neuron of the inhibitory population, 
scaled by a coefficient Wi e = 15, was the inhibitory input for the excitatory 
population. Excitatory/inhibitory connections are shown by solid/dashed lines. 
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these figures show the population activities computed by the EDMs and 
the bottom panels the activity derived from the direct simulation of 9,000 
IF neurons. The blue curves, scaled along the left axis, show the average 
spike rate per neuron in the population, the magenta and yellow curves, 
scaled along the right axis, show the excitatory and inhibitory external 
inputs to the population, respectively, and the magenta and yellow dashed 
curves, also scaled along the right axis, show the excitatory and inhibitory 
feedback inputs to the population. We should have performed the direct 
simulations with a larger number of IF neurons to obtain smoother spike 
rates and currents. Nevertheless, the activities computed by the EDMs 
are in close agreement with those derived from the direct simulation. 

2.3 Partial Conclusions 

Given a dynamical model of a neuron, we now know how to derive an EDM 
for a population of such neurons. For an IF model of a neuron, here we 
have shown that EDMs accurately approximate population activity (i.e., 
the pdf of the trans-membrane voltage, p(v,t), and the average firing 
rate per neuron, r(f)). The next step in this sub project is to estimate 
connectivity parameters (e.g., Wi e and W Kl in Figure 0 from simulated 
data. 
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Figure 8: The top and bottom panels show the activities of the excitatory 
population represented by the population equation, and by the simulation of 
9,000 IF neurons. The blue line, scaled on the left axis, plots the firing rate of the 
population. The magenta and yellow lines represent excitatory and inhibitory 
currents, respectively; and the solid and dashed lines represent external and 
feedback currents, respectively. The currents are scaled on the right axis. Note 
that the similarity between the average firing rate per neuron obtained from the 
population equation and from the direct simulation (blue lines in the top and 
bottom panels). 
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Figure 9: Same as Figure [8] but for the inhibitory population in the two- 
populations model. 
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3 Reducing dimensionality in EDMs 

In the previous section we showed that Ensemble Density Models (EDMs) 
accurately approximated the average firing rate per neuron and the proba¬ 
bility density function (pdf) of direct simulations of ensembles of integrate- 
and-fire (IF) neurons. In networks of EDMs, we want to estimate con¬ 
nectivity parameters and state variables (i.e., the pdfs of the different 
ensembles) from recorded ensemble firing rates. The state space of each 
previously reported EDM contained 210 variables. To make the estimation 
of parameters and state variables in networks of EDMs feasible/efficient it 
would be helpful to find low-dimensional approximations of EDMs. Here 
we report the approximation power of one such low-dimensional approxi¬ 
mation method. This method was inspired by the moving basis technique 
in iKnightl [200d ]. 

3.1 Method to find low-dimensional approxima¬ 
tions of EDMs 

The evolution of the ensemble pdf, p(v,t), is given by: 


p(v,t) = Q(s(t))p(v,t) 


(4) 


where Q(s(t)) is a differential operator that depends on the stimulus s(t). 
The normalized voltage v in Equation [I] ranges in the unit interval. To 
numerically solve this equation, we discretize v, {vt = i/N : 1 < i < N} 
and p, {pi(t) = p{vi — A/2, i) : 1 < i < IV}, giving a discretization of 
Equation U 


p{t) = Q(s(t))p(t) 


(5) 


where Q(s(t)) € R iVxJV is the matrix representation of the differential op¬ 
erator Q(s(t). Equation [5] is a system of N differential equations. The 
objective of the dimensionality reduction method described below is to 
approximate the evolution of p(t) using a system of M differential equa¬ 
tions, where M -C N. For this, the ensemble pdf p(t) is represented in a 
basis of eigenvectors of the differential matrix Q(s(t), where many coef¬ 
ficients in the new representation can be discarded without much loss in 
approximation power. 

3.1.1 Representing the ensemble pdf in a new basis 

Let {4>n(s(t)) : 0 < n < N} and (A n (s(t)) : 0 < n < N} be the eigenvec¬ 
tors and eigenvalues, respectively, of Q(s(t)): 


Q(s(t))<p n (s(t)) = A n(s(t))(f>n(s(t)) 


( 6 ) 


or in matrix notation: 


Q{s{t))$(s(t)) = $(s(i))A(s(t)) 


(7) 
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where 4> n {s(t)) is the nth column of the matrix <f>(s(f)) and A n (s(t)) is 
the nth diagonal element of the diagonal matrix A(s(t)). 

Assuming the eigenvectors are linearly independent, we represent p(t ) 
as: 


N 

P(t ) = ^2 a n (t)4>„(s(t)) (8) 

n= 1 


or in matrix notation: 


p(t) = $(s(f))a(t) 
From Equations [9] and |T] it follows: 


(9) 


Q(s(t))p{t) = Q(s(t))$(s(t)) a (t) = $(s{t))A(s(t)) a(f) (10) 


Using a backward difference to approximate the time derivative of p(t): 


= p(t) — p(t — At) 
; At 

in Equation [5] we obtain: 


( 11 ) 


P(t) - A tQ(s(t))p(t) = p(t - At) (12) 

Now applying Equations [9] and [TO] to Equation [12] we get: 


[$(«(/:)) (/ — AtA(s(t)))] a(t) = $(s(t — At)a(t — At) (13) 

or: 


a(t) = [(/ — AtA(s(t))) 1 $(s(t)) 1 $(s(t — At)] a(t — At) (14) 

Equation M provides the evolution of the coefficients a (t) in Equa¬ 
tion |2] It just expresses the evolution of the ensemble pdf in Equation 0 
in another basis. It is an exact formula (i.e., it is not an approximation) 
and has the same dimensionality as Equation [5] The importance of Equa¬ 
tion m is that one can approximate the evolution of the ensemble pdf by 
discarding many components of a(f), as we explain next. 

3.1.2 Reducing dimensionality in the new basis 

To understand why one can discard many coefficients in the representation 
of the ensemble pdf of Equation 0 without much loss in approximation 
power, consider the case where the stimulus, s(f), is constant. In such a 
case, the eigenvectors and eigenvalues will neither depend on the stimulus; 
i.e., $(s(t)) = <f> and A(s(f)) = A. Then Equations [5] and [5] reduce to: 
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and: 


Pit) = Qp(t) 


( 15 ) 


pit) = <f> a [t] 

Taking derivatives in Equation [TB] we obtain: 


(16) 


p(f) = $a(t) (17) 

and substituing Equation [16] in Equation [T5] and applying Equation 0 we 
get: 


p(i) = Q$a(t)=Ma(i) (18) 

Equating the right hand sides in Equations 1171 and 1181 and pre multi¬ 
plying by <4> _1 , gives: 


a(() = Aa(t) (19) 

or: 

di(t) — Xi cn(t) (20) 

with solution: 

ciiit) = exp(A it) aii0 ) (21) 

Thus, the absolute value of a low dimensional coefficient at evolves as: 


|a-£(£)| = exp(K(Ai)t) |ai(0)| (22) 

Because all eigenvalues have negative real part (with the exception of 
the zero eigenvalue), the coefficients associated with non-zero eigenvalues 
will decay to zero, and the speed of this decay will be proportional to the 
absolute value of the real part of the corresponding eigenvalue. Therefore, 
to achieve dimensionality reduction in Equation [14] we may discard those 
coefficients associated with eigenvalues with larger absolute value of their 
real part, since these coefficients will rapidly decay to zero. 
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3.2 Evaluation of the method 


We first study how well low-dimensional EDMs approximate the average 
firing rate per neuron in direct stimulations (Section 13.2. Il l and then how 
they approximate the ensemble pdf (Section l3.2.211 . 

These studies were peformed in data simulated from the network of 
EDMs illustrated in Figure [3 The network was driven by an excitatory 
sinusoidal input to the excitatory population (shown by the dotted curve 
and scaled along the right axis in the top panel of Figure [Toll . The average 
firing per neuron in this population, scaled by a constant W Kt = 50, drove 
the inhibitory population (shown by the dotted curve and scaled along the 
right axis in the top panel of Figure fTIll . In turn, the average firing rate 
per neuron in the inhibitory population, scaled by a constant Wi e = 15, 
inhibited the excitatory population. Both populations contained feedback 
(i.e., each neuron received 10 inputs from ten other cells in the same 
population and 80% of these inputs were inhibitory). 

3.2.1 Firing Rates 

The top panels of Figures [10] and [11] show the average firing rate per neu¬ 
ron obtained from direct simulation (grey curve), from a full-dimensional 
EDM (red curve), and from low-dimensional EDMs (the blue, cyan, and 
red curves correspond to 17, 5, and 1 moving basis, respectively). The 
full dimensional EDM and its low dimensinal approximations with 17 and 
5 moving basis almost perfectly approximate the average firing rate per 
neuron of the direct simulation. The approximation power of the EDM 
with only one moving basis is not as good, but it look reasonable. 

3.2.2 Ensemble Proability Density Functions 

We compared the normalized histogram of number of directly simulated 
neuron per voltage bin (i.e., the ensemble pdf from direct simulation) with 
the pdfs calculated with EDMs (i.e., Equation[5]). For this we computed at 
every time step the Kullback-Leibler (KL) divergence (in bits) between the 
pdfs obtained by direct stimulation and those obtained from EDMs. These 
KL divergences are shown in the bottom panels of Figures flOl and [111 

We see that the pdfs obtained from EDMs were good approximation of 
the pdfs from the direct simulation at times of large average firing rate per 
neuron (between 0.73 and 0.85 seconds and between 1.03 and 1.2 seconds 
in the bottom panel of Figure flOl and between 0.72 and 0.82 seconds and 
between 1.03 and 1.18 seconds in the bottom panel of Figure fTIll . 

At times of low averaged firing rate per neuron the difference between 
the pdfs obtained from direct simulation and those obtained from EDMs 
were an order of magnitude larger for the inhibitory than for the excitatory 
ensemble. This is probably because the excitatory ensemble was driven 
by a large and smooth sinusoidal input while the inhibitory ensemble 
was driven by the weaker and non-smooth average firing rate per neuron 
of the excitatory ensemble. Also, at most time points, the larger the 
number of moving basis in low-dimensional EDMs, the better the EDM 
pdf approximated the pdf obtained from direct simulation, as can be seen 
more clearly in the bottom panel of Figure [IT] 
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Figure 10: Average firing rate per neuron (top) and KL divergence between 
the ensemble pdf calculated by direct simulation and that caculated by EDMs 
(bottom) for the excitatory ensemble. 
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Figure 11: Average firing rate per neuron (top) and KL divergence between 
the ensemble pdf calculated by direct simulation and that caculated by EDMs 
(bottom) for the inhibitory ensemble. 
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Figure 12: Ensemble pdfs for the inhibitory ensemble at 0.81s. 


To try to understand why the pdfs obtained from direct simulation 
were different from those obtained from low-dimensional EDMs, we plot¬ 
ted these pdfs for the inhibitory ensemble in an interval of low average 
firing rate per neuron (from 0.81 seconds in Figure fl2l to 1.09 seconds in 
Figure l26l) . In the transition between weak to zero average firing rate 
per neuron (from 0.81 seconds in Figure IT2l to 0.93 seconds in Figure fl8l) 
the low-dimensional pdfs moved faster towards lowers voltage than the 
pdfs from the full-dimensional EDM and those from direct simulation. 
Similarly, in the transition between zero to weak average firing rate per 
neuron (from 0.95 seconds in Figure ITiTI to 1.09 seconds in Figure E51l the 
low-dimensional pdfs moved faster towards higher voltages than the pdfs 
from the full-dimensional EDM and those from direct simulation. This 
suggests that the moving basis discarded in the low-dimensional approxi¬ 
mations of EDMs help prevent the EDM pdf to transition too fast to and 
away from the pdf corresponding to large average firing rate per neuron. 

3.3 Partial conclusions 

We have developed a method to reduce the dimensionality of the state 
space of EDMs. We observed that the pdfs from low-dimensional EDMs 
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Figure 13: Ensemble pdfs for the inhibitory ensemble at 0.83s. 
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Figure 14: Ensemble pdfs for the inhibitory ensemble 0.85s. 
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Figure 15: Ensemble pdfs for the inhibitory ensemble 0.87s. 
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Figure 16: Ensemble pdfs for the inhibitory ensemble 0.89s. 
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Figure 17: Ensemble pdfs for the inhibitory ensemble 0.91s. 
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Figure 18: Ensemble pdfs for the inhibitory ensemble 0.93s. 
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Figure 19: Ensemble pdfs for the inhibitory ensemble 0.95s. 
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Figure 20: Ensemble pdfs for the inhibitory ensemble 0.97s. 
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Figure 21: Ensemble pdfs for the inhibitory ensemble 0.99s. 
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Figure 22: Ensemble pdfs for the inhibitory ensemble 1.01s. 
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Figure 23: Ensemble pdfs for the inhibitory ensemble 1.03s. 
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Figure 24: Ensemble pdfs for the inhibitory ensemble 1.05s. 
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Figure 25: Ensemble pdfs for the inhibitory ensemble 1.07s. 
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Figure 26: Ensemble pdfs for the inhibitory ensemble 1.09s. 
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move faster to and away from the high-average-firing-rate-per-neuron pdf 
than the pdf derived from direct simulation or that from full EDMS. 
However, these differences did not have a large impact on the average 
firing rate per neuron produced by low-dimensional EDMs, which almost 
perfectly approximates the average firing rate per neuron computed by 
direct simulation. 

Our next step is to estimate connectivity parameters and state space 
variables in networks of EDMs from recorded spike rates. We will first do 
this with simulated data. 
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G=5, f=l 


G=5, f=0 


Figure 27: Network of EDMs. 


4 Estimating parameters of EDMs 

Here we evaluate a maximum likelihood method to estimate the connec¬ 
tivity parameters 9 = [W e i,Wi e ] in the EDMs network in Figure [27] We 
seek the connectivity parameters for which a set of N pairs measurements, 
Yn = [y(0),. .., y (N — 1)], are most probable; ie: 


9 m i = arg max P(Y N \0) 
e 

The pair of measurements at time n, y (n) comprises measurements from 
the excitatory and inhibitory ensembles; i.e., y(n) = [j/ e (n), j/i(n)]. 

4.1 Noise Model 

As a first approximation we make the following three assumptions on the 
noise of the measurements . These are the assu mptions made in previous 
work by our collaborators [Kostuk et all 12012 ]. 

1. Noise is independent in time; i.e., P(Yn\0 ) = n„ =1 P(y(n)|0). 

2. Noise is independent in the different populations; i.e., P(y(n)|0) = 
P{y e {n)\9)p[yi(n)\B). 

3. Gaussian noise with a known variance, <r 2 , in each population; i.e., 
P(y.(n)\6) = N(y.(n)\r.(n\9), a 2 ), where r.(n\6) is the activity gen¬ 
erated by the ensemble at time n. 

With these assumptions the log likelihood function of the model pa¬ 
rameters reduces to: 


logP(y JV |0) = A-E^ =l 


(Ve(n) 


r e (n\6)) 2 + ( y e (n ) 
2a 2 


r e (n\6)) 2 


The red lines in Figure [28] plots the activity generated by the excita¬ 
tory and inhibitory ensemble with connectivity parameters W e i = 50 and 
Wie = 15. The blue lines plot the noisy measurements (a = 2) that we 
use below for parameter estimation. 
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Figure 28: EDMs spike rates (red) and noisy measurements (blue, a = 2) used 
to estimate the connectivity parameters in the network of Figure [571 
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4.2 Optimization Surface 

The levelplot in Figure [29] shows the optimization surface for the set of 
parameters show in the axes. We see that it has a convex shape with peak 
at the true parameters. Thus, an iterative gradient-ascent optimization 
procedure should climb to the maximum-likelihood parameters from any 
starting set of parameters. 

4.3 Gradient of Log-Likelihood Function 

To compute the gradient of the log-likelihood function at a given set of 
parameters 9 = [w e i, Wi e ] one needs to integrate the EDMs to generate 
at each time step the activities of both ensembles, r e (n,9) and ri(n,9), 
as well as the ensembles pdfs, p e (n\0) and pi(n\6). With these quanti¬ 
ties in hand, the gradient of the log-likelihood function can be computed 
recursively in a second integration step fEquation 1231) . 


Slog P(Y m \0) 

8Wei 
dr e [n; 9} 
dW ei 
dp e [n\9] 
dWei 
dQ e [n ; 9] 
dWei 


E ra 
T?, = l 


{y e [n} — r e [n: 9]) dr ^ e] + ( yi[n\ - n[n; 0]) dr ^ e] 


A E rw dp e [n\9] 

Ava * [n] ^^w— ] 

dQ e \n — 1; 01 . . , . r dpAn — 110] 

At Ve L.. l Pe[n - 1| 9} + [/ + A tQ e [n - 1; 0]] -LT 


= -W e 


8We; 

dn[n-l\9] (2) 

1 dWei 


OWe; 


(23) 


The white arrows in Figure [29] point in the direction of the gradi¬ 
ent, and are scaled according to the gradient magnitude. Note that, as 
expected, the gradient direction is perpendicular to the level lines. 

With a convex log-likelihood function, for which we can compute its 
gradient, we can now use an iterative gradient ascent procedure to max¬ 
imize this function. The green line in Figure [2H shows a gradient-ascent 
trajectory that in only three steps accurately approximated the maximum 
likelihood parameters. 
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Figure 29: Log-likelihood function, its gradient, and a sample gradient ascent 
path. The contour plot plots the log-likelihood function for the parameter values 
shown on the axes for the network in Figure 1271 The white arrows show the 
log-likelihood gradient computed analytically. The green curve is a gradient 
ascent trajectory starting at W e i = 37 and Wi e = 3. 
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5 Modeling evoked auditory activity in 
rodents with EDMs 


Here we show that the average firing-rate per neuron predicted by EDMs 
well approximates the evoked high-gamma power (HGP) from auditory 
neurons in response to stimulation with pure-tones. 

5.1 Recordings 

We used very high-resolution simultaneous surface and laminar recordings 
(Figure I3U11 in anesthetized rodents stimulated with pure tones at different 
frequencies and amplitudes. The blue trace in Figure I3T1 shows the high- 
gamma power (HGP) evoked by the presentation of tones in one sample 
electrode of the array. Each vertical dotted color line marks the onset of 
a tone of a corresponding frequency. 

We see that these evoked waveforms are very stereotypical. Some 
waveforms have a large peak followed by a bump. Other waveforms have 
two peaks, with a smaller peak preceding a larger one. We also see inf 
Figure [32] waveforms that have more than one peak preceding a larger one 
and waveforms having a larger peak followed by a smaller one. 

5.2 Qualitative analysis 

The recorded HGP waveforms appear remarkably similar to those pro¬ 
duced by EDMs when stimulated by sinusoids. To confirm this similar¬ 
ity we manually chose parameters for EDMs to reproduce the shape of 
recorded waveforms (Figures 13311341 1351 and ED- 

5.3 Quantitative analysis 

Using a similar method as described above to learn connectivity param¬ 
eters in networks of EDMs, we etimated parameters so that an EDM 
approximated as close as possible the HGP evoked by the presentation of 
pure tones to rodents. In total we estimated 13 parameters: three param¬ 
eters of an EDM, two parameters for its initial condition, four parameters 
for the sinusoidal excitatory input, and four parameters for the inhibitory 
input. The recorded and approximated HGP are given by the blue and 
red solid curves, respectively, in Figure [37] 

5.4 Partial Conclusions 

These preliminary results show that EDMs, besides accurately approxi¬ 
mating the average firing rate of ensembles of simulated IF neurons, well 
approximate the HGP in the auditory cortex of rats evoked by auditory 
tones. 
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Figure 30: Electrodes used for simultaneous surface and laminar recordings. 
Surface recordings were obtained from an 8 x 8 grid of subdurally implanted 
electrodes covering an area of 1.6 mm 2 . Laminar recordings were obtained from 
a 32-channel politrodes of length 650 /*m. These rercordings were combined 
with optogenetic manipulations. 
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Figure 31: High-gamma power evoked by the presentation of pure tones at 
different frequencies and amplitudes. Blue traces show the HGP evoked by the 
presentation of tones in one sample electrode of the array. Each vertical dotted 
color line marks the onset of a tone of a corresponding frequency. 
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Figure 32: More examples of HGP evoked by the presentation of pure tones. 
The format is as in Figure ED 
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Figure 33: The blue curve in the inset shows the average firing-rate simulated 
by an EDM with manually chosen parameters. The magenta and yellow solid 
curves plot its excitatory and inhibitory inputs, respectively. The first waveform 
is qualitatively similar to the ECoG waveforms with a bump following a large 
peak. 
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Figure 34: The blue curve in the inset shows the average firing-rate simulated by 
an EDM with manually chosen parameters. The magenta and yellow solid curves 
plot its excitatory and inhibitory inputs, respectively. The second waveform is 
qualitatively similar to the ECoG waveforms with a lower peak preceding a large 
one. 
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Figure 35: The blue curve in the inset shows the average firing-rate simulated 
by an EDM with manually chosen parameters. The magenta and yellow solid 
curves plot its excitatory and inhibitory inputs, respectively. EDMs can generate 
waveforms with multiple lower peaks preceding a larger one. 
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Figure 36: The blue curve in the inset shows the average firing-rate simulated 
by an EDM with manually chosen parameters. The magenta and yellow solid 
curves plot its excitatory and inhibitory inputs, respectively. EDMs can generate 
waveforms with a larger peak preceding a smaller one. 
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[ 10.10648747 21.14308782 0.80066418 12.04027179 11.95724926 

722.23861385 279.88373368 4.54564969 - 1.16267511 634.38337831 



Figure 37: Learning ensemble properties to approximate physiological record¬ 
ings. We estimated 13 parameters so that the average firing rate per neuron 
predicted by an EDM approximates as close as possible the recorded HGP in 
response of a pure tone. The title shows the values of these parameters. The 
blue and red curves in the top panel plot the recorded HGP and the predicted 
average firing rate per neuron, respectively. The dotted green and magenta 
curves in the top panel shows the estimated excitatory and inhibitory inputs to 
the EDM. The curve in the bottom panel plot the estimated initial condition 
probability density function of the EDM. The average firing rate predicted by 
the EDM is a good approximation of the recorded HGP. 
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6 Conclusions 


We have shown that an ensemble model accurately reproduce the probabil¬ 
ity density function of the transmembrane voltage, as well as the average¬ 
firing rate per neurona, in a large ensemble of integrate-and-hre simulated 
neurons (Section [3J. We developed and evaluated methods to reduce the 
dimensionality (Section [3J) and estimate parameters (Section [0. Finally 
we demonstrated the feasibility of EDMs to model the high-gamma power 
evoked by pure tones in the auditory cortex of rodents. 

The possiblity of quantitatively model the activity of ensemble of neu¬ 
rons may allow us to uncover fundamental computations performed by 
neural ensembles. 
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